Ground state of the Bethe-lattice spin glass and running time of an exact 

optimization algorithm 
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We study the Ising spin glass on random graphs with fixed connectivity z and with a Gaussian 
distribution of the couplings, with mean jj, and unit variance. We compute exact ground states 
by using a sophisticated branch-and-cut method for z — 4,6 and system sizes up to 1280 spins, for 
different values of /j.. We locate the spin-glass/ferromagnet phase transition sd, ^ — 0.77±0.02 {z — 4) 
and fj, = 0.56 ± 0.02 {z = 6). We also compute the energy and magnetization in the Bethe-Peierls 
approximation with a stochastic method, and estimate the magnitude of replica symmetry breaking 
corrections. Near the phase transition, we observe a sharp change of the median running time of 
our implementation of the algorithm, consistent with a change from a polynomial dependence on 
the system size, deep in the ferromagnetic phase, to slower than polynomial in the spin-glass phase. 



I. INTRODUCTION 

Recent years have seen an increasing interaction be- 
tween the fields of combinatorial optimization and sta- 
tistical physicsiiSi^. On one hand, several problems in 
the statistical physics of disordered systems have been 
mapped onto combinatorial problems, for which fast com- 
binatorial optimization algorithms are available"*'^. This 
has provided valuable insights into questions that are 
hard to investigate with traditional techniques, such as 
Monte Carlo simulations. On the other hand, concepts 
and methods from statistical physics are increasingly ap- 
plied to combinatorial optimization^. 

Easy/hard thresholds analogous to phase transitions 
have been observed in random instances of optimization 
and decision problems, including satisfiability (SAT)^-'^ , 
vertex-cover^ {VC), number partitioning^., and others. 
There is currently much interest in understanding how 
phase transitions affect the performance of combina- 
torial algorithms, following the observation^^ that the 
averagei^ or typical (i.e. median) running time of some 
algorithms exhibits a sharp change in the vicinity of a 
phase transition. For example, in 3SAT, a special case 
of SAT, and in VC, the typical running time of exact 
backtracking algorithms change o^^d^ from a polynomial 
dependence on the input size in the "solvable" region, to 
exponential dependence in the " unsolvable" region. This 
provides an insight into the performance of algorithms 
that goes beyond the worst-case running time tradition- 
ally considered in complexity theory. (Note, however, 
that from the behavior of individual algorithms, strictly 
speaking, one cannot draw conclusions about the "typ- 
ical hardness" of a problem itself). Recently, statistical 
physics techniques have been fruitfully applied to study 
easy/hard transitions and algorithmic performance^. 

In this paper, we apply a branch-and-cut algorithm, 
a technique developed in combinatorial optimization, to 
find the ground state of the Ising spin glass on random 
graphs with fixed coordination number (also called Bethe 



lattices) . 

Our motivation is twofold. The first goal is algo- 
rithmic: we want to characterize the typical running 
time of our algorithm, notably its behavior across the 
zero-temperature spin-glass/ferromagnet phase transi- 
tion that occurs when varying the mean of the random 
couplings. The interest of this stems from the importance 
of branch-and-cut as a general technique in combinatorial 
optimization, and from the fact that finding the ground 
state of a spin glass is a prominent example of a hard 
optimization problem arising from statistical physics (in 
general, it is NP-hardU, see Section IIII|I . The perfor- 
mance of branch-and-cut for this application has not been 
investigated in detail before (sec, however, Refs. 0, IT^ 
and 17), and here we fill this gap. An interesting aspect is 
that, unlike in SAT, VC and other classical combinatorial 
problems, here averaging over random instances is physi- 
cally motivated. To our knowledge, the only other study 
relating a "physical" phase transition to algorithmic per- 
formance is that of Middletoni^, which investigates the 
typical running time of the matching algorithm for the 
random-field Ising model, which however is polynomial 
everywhere in the parameter space. 

We find that the median running time of our algorithm 
varies sharply near the spin-glass/ferromagnet transition, 
indicating a change from polynomial time deep in the fer- 
romagnetic phase, to slower than polynomial in the spin 
glass phase. We also observed a similar behavior for spin 
glasses on regular lattices in two and three dimensions, 
but will not report it here. 

The second motivation for the present work lies 
in the ground-state properties of the Bethe-lattice 
spin glass, which recently have attracted a renewed 
interesti2i22i2ii2^. Using branch-and-cut, we compute the 
ground state energy and magnetization, and locate the 
spin-glass/ferromagnet phase transition. This provides 
a useful test of recently developed analytical methods 
to treat diluted spin glass modelsi^iSLS^. We solve the 
model in the Bethe-Peierls (BP) approximation (equiv- 
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alent to the replica symmetric approximation in the 
replica formalism) using a stochastic approach proposed 
m By comparing the branch-and-cut re- 

sults with the BP results, we estimate the magnitude of 
the replica symmetry breaking corrections to the ground 
state energy and magnetization, finding that they are 
small. 

The rest of the paper is organized as follows. In Sec- 
tional we introduce the Bethe-lattice spin glass model. 
In Section Hm we describe the branch-and-cut algorithm 
used to calculate the exact ground states of the model. In 
Section Hvl we describe the Bethe-Peierls approximation 
and the stochastic procedure used to solve it. In Sec- 
tion we present our branch-and-cut and BP results 
for the ground state energy and the zero temperature 
phase transition. In Section IVII we show that this tran- 
sition coincides with a change of the typical running time. 
Finally, Section [VIII summarizes our results. 

II. MODEL 

The system considered consists of N Ising spins Si = 
±1 sitting on the nodes of a graph G = {V,E), where 

V — {1, . . . , N} is the set of nodes and E — {(«, j)} C 

V xV is the set of edges of the graph. The energy of the 
system is given by 

where the couplings Jij are independent, identically dis- 
tributed random variables drawn from a Gaussian distri- 
bution P(-) with mean /i and unit variance, 

P(J) = ^exph(J-/.)V2] (2) 

We consider the case in which G is a random graph 
with fixed connectivity z, or z-regular graph, where each 
spin interacts with exactly z neighbors. This provides 
a convenient realization of a Bethe lattice, which avoids 
some complications associated to the usual construction 
of a Bethe lattice from a Cayley treei^. Frustration is 
induced by large loops, the typical size of a loop being 
of order log(A^). Small loops are rare, giving the graph a 
local tree-like structure, and therefore the mean field ap- 
proximation is exact. A related model is the Viana-Bray 
model2^, in which the connectivity is a Poisson variable 
with finite mean. Finite-connectivity or "diluted" mod- 
els provide a better approximation to finite-dimensional 
spin glasses than the infinitely-connected Sherrington- 
Kirkpatrick model. Furthermore, they are directly re- 
lated to optimization problems such as graph partition- 
ing and coloring. 

Although it has long been known that replica sym- 
metry is broken in these two modelsSS*2Si2i, until re- 
cently a replica symmetry broken solution could be 
found only in some limit cases. Mezard and Parisi 



recently introducedi2i2£ a "population dynamics" algo- 
rithm which allows a full numerical solution at the level 
of one step of replica symmetry breaking. Explicit re- 
sults were derived for the Bethe-latticc spin glass with 
the symmetric ± J disorder distribution, but not for the 
Gaussian distribution considered here or for a non-zero 
mean. Previous numerical studies of this model can be 
found in Refs. 22 28 29. For a complete discussion of 
the Bethe-lattice spin glass, see Ref. 19: and references 
therein. 



III. BRANCH-AND-CUT ALGORITHM 

The problem of finding a ground state of the Hamil- 
tonian in Eq. is in general computationally demand- 
ing. For a generic graph G, it is NP-Ziarrf'^'^'*. For NP- 
hard problems, currently only algorithms are available, 
for which the running time increases faster than any poly- 
nomial in the system size, in the worst case (see Section 
VI for a brief description of complexity classes). In the 
special case of a planar system without magnetic field, 
e.g. a square lattice with periodic boundary conditions in 
at most one direction, efficient polynomial-time matching 
algorithms'^" exist. For the square lattice with periodic 
boundaries in both directions, polynomial algorithms ex- 
ist for computing the complete partition function for the 
± J distribution"^^ and for the case in which the coupling 
strenghts are bounded by a polynomial in the system 
size^^. In practice, both algorithms can only reach rela- 
tively small system sizes. 

For the Bethe lattice considered here (and for reg- 
ular lattices in dimension higher than two), no poly- 
nomial algorithm is known. Heuristic algorithms 
recently used include simulated annealing'^-, "multi- 
canonical" simulations!, genetic algorithma^^*^, ex- 
tremal optimization^^, a hierarchical renormalization- 
group based approach^!, and the cluster-exact approx- 
imation algorithm^. Here, however, we are interested 
in investigating the running time of an exact, determin- 
istic algorithm, since in this case the running time to 
find the exact ground state is a well defined quantity. 
We study the branch-and-cut methodi^i^S (see Ref. 113 
for a tutorial on optimization problems and techniques, 
including branch-and-bound and branch-and-cut), which 
is currently the fastest exact algorithm for computing 
spin glass ground statesii, with the exception of the 
polynomial-time special cases mentioned above. As the 
branch-and-cut method is basically branch-and-bound 
with cutting planes, we also did some experiments with 
a pure branch-and-bound algorithm^ii^, which however 
can only deal with much smaller system sizes, finding a 
qualitatively similar, but less pronounced, variation of 
the running time across the transition. 

In the rest of this Section, we repeat a short description 
of the branch-and-cut method already given in Ref.m 
to the benefit of the reader. It is convenient to map the 
problem of minimizing the Hamiltonian in Eq.(^| into a 
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maximum cut problem. Consider a graph G = [V^E)^ 
and let assign weights {Kij} to the edges. Given a parti- 
tion of the node set V into a subset W CV and its com- 
plement V \ W, the cut S{W) associated to W is the set 
of edges with one endpoint in W and the other endpoint 
mV\W, namely S{W) = {{ij) e E \ i e W,j e V\W}. 
The weight of S{W) is defined as the sum of the weights 
of the cut edges, J2{ij)e5iw) -^u- maximum cut is 

a node partition with maximum weight among all parti- 
tions. It can be shown'^^ that minimizing the Hamilto- 
nian in Eq.(^ is equivalent to finding a maximum cut of 
G with the assignment Kij = —Jij- 

The branch-and-cut algorithm solves the maximum cut 
problem through simultaneous lower and upper bound 
computations. By definition, the weight of any cut gives a 
lower hound on the optimal cut value. Thus, we can start 
from any cut and iteratively improve the lower bound us- 
ing deterministic heuristic rules (local search and other 
specialized heuristics, see Ref. |43| for details). How do 
we decide when a cut is optimal? This can be done by 
additionally maintaining upper hounds on the value of 
the maximum cut. Upon iteration of the algorithm, pro- 
gressively tighter bounds are found, until optimality is 
reached. 

Since the availability of upper bounds marks the dif- 
ference between a heuristic and an exact solution, we 
now summarize how the upper bound is computed (for 
more details, see Ref.lisll To each edge {ij) we associate 
a real variable Xij and to each cut 5{W) an incidence 
vector -x^^^^f e IR-^ with components x^j^^ associated 
to each edge {ij), where x\f^^ = 1 if {ij) € 6{W) and 

Xij^^ = otherwise. Denoting by Pc{G) the convex 
hull of the incidence vectors, it can be shown that a hasic 
optimum solution^ of the linear program 

max{ ^ JijXij I X e Pc{G)}. (3) 

is a maximum cut. In order to solve with linear pro- 
gramming techniques we would have to express Pc{G) in 
the form 

Pc{G) = {x e -R^ \ Ax < h,0 < X < 1} (4) 

for some matrix A and some vector b. Whereas the exis- 
tence of A and h are theoretically guaranteed, even sub- 
sets of Ax < h known in the literature contain a huge 
number of inequalities that render a direct solution of 
(PJ impractical. 

Instead, the branch-and-cut algorithm proceeds by op- 
timizing over a superset P containing Pc (G) , and by it- 
eratively tightening P, generating in this way progres- 
sively better upper bounds. The supersets P are gener- 
ated by a cutting plane approach. Starting with some P, 
we solve the linear program maxj^^^^^^^ JijXij \ x G P} 
by Dantzig's simplex algorithmic. Optimality is proven 
if either of two conditions is satisfied: (i) the optimal 



value equals the lower bound; (ii) the solution vector x 
is the incidence vector of a cut. 

If neither is satisfied, we have to tighten P by solving 
the separation prohlem. This consists in identifying in- 
equalities that are valid for all points in Pc{G), yet are 
violated by x, or reporting that no such inequality ex- 
ists. The inequalities found in this way are added to the 
linear programming formulation, obtaining a new tighter 
partial system P' C P which does not contain x. The 
procedure is then repeated on P' and so on. 

At some point, it may happen that (i) and (ii) arc not 
satisfied, yet the separation routines do not find any new 
cutting plane. In this case, we hranch on some fractional 
edge variable Xij (i.e. a variable Xij ^ {0,1}), creating 
two subproblems in which Xij is set to and 1, respec- 
tively. We then we apply the cutting plane algorithm 
recursively for both subproblems. 

IV. BETHE-PEIERLS APPROXIMATION 

We recall here the zero-temperature formulation of the 
BP approximation, loosely following Ref.[23. We consider 
the Hamiltonian Eq.(^ on a random graph with fixed 
connectivity z = fc + 1, in which however some spins Si 
{cavity spins) have only k neighbors. The random cou- 
plings are drawn from a distribution P{J). The BP ap- 
proximation consists in assuming that the ground state 
energy of this system is given by £^ = const. — hiSi, 
where the sum runs over all cavity spins. The cavity fields 
hi, implicitly defined by this relation, are independent, 
identically distributed random variables when considered 
as a function of the random couplings. Their distribution 
P{h) is the central object of interest, and satisfies a re- 
cursion relation derived as follows. Suppose we add a new 
spin 5*0 to the system, which interacts with k pre-existing 
cavity spins Si, . . . ,Sk through couplings Ji, . . . , Jk, and 
we minimize the energy with respect to Si, . . . , Sk- Now 
Sq is a cavity spin, and it is easily shown that its cavity 
field ho is given by 

k 

ho = ^u{Ji, hi) , (5) 

i=l 

where u{Ji, hi) — ^ {\hi + Ji\ — \hi — Jj|). This provides 
a recursion relation for P{h) as the J's fluctuate accord- 
ing to P(J). 

Given a spin Sq interacting with k + 1 neighbors with 
couplings Ji, . . . , Jk+i, the internal field H acting on Sq 
is 

k+l 

H = Y,<'h,h.d. (6) 

i=l 

Therefore, if we know P{h) we can determine the proba- 
bility distribution P{H), and the magnetization 

mi3P = lim I dH P{H)sgn{H) , (7) 
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where sgn(x) is the sign function and e is a smaU field that 
breaks the symmetry of Eq. Q with respect to changing 
the sign of all cavity fields. 

The knowledge of P{h) is also sufficient to determine 
the ground state energy of the system. As shown in 
Ref. 123, this can be expressed as 

esp = [AE('^-^[AE('^], (8) 

where [• • •] is the expectation value with respect to P{ J) 
and P{h), and the quantities AE^^\ AE'^^^ are given by 

fc+i fe+i 
Ai;(i) =-^a(J„/iO- |^^^(J»>»)| (9) 

2=1 i=l 

and 

A£;(2) = - max(/ijS'j + hjSj + J^jS^SJ) , (10) 

where a{Ji,hi) = ^ {\hi + J^] + \hi - J,\) and Si,Sj in 
Ea. (|10|l are two randomly chosen cavity spins which we 
connect with the coupling Jy . 

The BP recursion, especially at finite temperature, 
has been studied extensively (see Ref. and references 
therein). In particular, Mezard and Paris:^° have given 
an analytic expression of P{h) for a binary P{J). Klein 
et ali^ solved the finite temperature BP recursion for 
Gaussian couplings with /i = in the vicinity of the 
spin-glass/paramagnet transition. No analytical solution 
has been derived for Gaussian couplings at T = 0, to 
our knowledge, although Klein et al,*^ derived an ana- 
lytic solution near the spin-glass/ferromagnet transition 
/i = /ic within the mean random-field approximation. 

Here, we employ the stochastic iterative procedure pro- 
posed by Mezard and ParisiSS for the more general one- 
step replica symmetry broken case (see also Ref. for a 
previous application of a similar method). We consider 
a population of Af sites, to which we associate Af cav- 
ity fields, which are initially assigned at random (with a 
small positive bias). We then select k sites at random, 
extract k couplings from P{J), compute ho from Eq.(|SJ), 
and assign ho as the new cavity field of a randomly chosen 
site. We iterate this procedure A4 times per site. After 
a certain number of iterations, the distribution P{h) will 
satisfy Eq.©. At each iteration, by merging A: -I- 1 ran- 
domly chosen sites we compute the internal field H with 
Eq.®, and AE'^^'> with Eq.l^, and by merging two sites 
we compute A£:(2) with Eq.lUnj. After discarding the 
first iterations, by averaging sgn(_ff), Ai?^^) and 

AE^^") over the remaining iterations we compute the es- 
timates oivTLBP and epp from Eqs.Q and (jSJ, and their 
statistical error from a binning procedure. We repeated 
the procedure for many values of /i, choosing AA — 10'' 
and M between 10^ and 10^, the larger population being 
for jjL near the transition point /ic- With TV = 10^, the 
iteration requires about one hour of computer time. 

We note that the BP approximation is known to be 
wrong, being equivalent to the replica symmetric solu- 
tion which is unstable. We have not attempted to use 



the generalization of the above procedure to one step of 
replica symmetry breakingiS, since for the Gaussian case 
considered it would require a significant computing time, 
and since the BP approximation gives sufficiently accu- 
rate results for our purposes. 



V. RESULTS 

We have studied the Ising spin glass on random graphs 
with fixed connectivity z = 4 and z = 6. The instance 
generator first builds a random regular graph with the 
algorithm described in Ref. We then assign the cou- 
plings Jij according to the distribution P{J) in Eq.Q. 

Using the branch-and-cut approach we were able to 
study graph sizes up to = 400 for z = 4 and /i < 0.9, 
and up to A^ = 200 for z — Q and /i < 0.7. For larger 
values of /i, we considered sizes up to A^ = 1280. The 
ground states for the smallest systems can be obtained 
within a second, while the longest computations lasted at 
most one day on a typical workstation. Incidentally, for 
Ising spin glasses on a regular grid, specialized heuristics 
exist that exploit the grid structure, making it possible 
to consider larger system sizes than for the model re- 
ported here. More one timing issues, in relation to the 
spin-glass/ferromagnet phase transition, is presented in 
Section Ell 

All the results were averaged over many samples (a 
sample, or instance, is a realization of the random graph 
with a realization of the couplings). The largest number 
of samples were considered in the vicinity of the phase 
transition, where the fluctuations of the magnetization 
are larger. Near the transition, for sizes N < 240 (z = 4) 
and A^ < 160 {z = 6) we computed around 5000 samples 
for each value of fi; for A^ = 400 (z = 4) and A^ = 200 
(z = 6), around 500 samples for each value of /i. For 
sizes larger than these, we computed up to 280 samples 
for each /i. In the following analysis of the ground state 
energy and magnetization, we consider only sizes up to 
A^ = 400 (z = 4) and A^ = 200 (z = 6), since for larger 
sizes the statistical error is quite large. In the analysis of 
running times we will include all sizes. 



A. Ground state energy 

We start by showing, in Fig. ^ the average ground 
state energy E{^,N), divided by zA^, as a function of 
fi for z = 4, 6 and two different system sizes. For suf- 
ficiently large /i, the system is completely magnetized, 
therefore the ground state energy depends linearly on fi, 
E{iJL, N) /N ~ z/i, as visible in the figure. For small /i the 
system is frustrated, hence the energy saturates. Note 
that E{0, N) scales as y/z, not as z, therefore the two 
curves diverge at small ^. The lines in Fig. ^ represent 
the numerical solution of the BP recursion obtained with 
a population size Af = 10'^ (we verified that with Af = 10^ 
the results are unchanged) and A4 = 10* iterations of the 
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FIG. 1: Normalized average ground state energy as a func- 
tion of the mean coupling strength, /i. The symbols represent 
the results of the branch-and-cut algorithm. Their statistical 
errors are smaller than the symbol sizes. The lines represent 
the results of the BP recursion. They are obtained by con- 
necting points spaced by A/i = 0.005 (A/i — 0.001 near the 
transition). Their statistical error is comparable to the line 
thickness. 

stochastic algorithm. Clearly, the branch-and-cut results 
agree well with the BP approximation. 

We extrapolate the branch-and-cut results to TV = oo 
by fitting the data with the form E/N ^ eoo + bN''^/^. 
As shown in Fig. |21 the finite size corrections are well de- 
scribed by a iV^^/'^ dependence for small /j,, although an 
A''"'^ correction fits reasonably well the data for other val- 
ues of oj between 0.6 and 1 as well. For large /z, the finite 
size corrections are very small. A N~^^^ correction was 
also found to fit well the numerical data by Boettcher^^, 
who computed the average ground state energy of the 
±J model for z up to z = 26 and N up to A = 2048 
using a heuristic algorithm. In Ref. lltA the finite-size 
dependence of the energy at T = 0.8 was studied, for 
the ± J distribution and z = 6, finding a finite-size expo- 
nent u! = 0.767(8), not far from 2/3. For the Viana-Bray 
model with fluctuating connectivity with mean z = 6, the 
value u! = 0.62 ±0.05, compatible with 2/3, was found**^, 
also using a heuristic algorithm. 

Fig. 1^1 also shows that the extrapolated energy, eoo, is 
very close to the BP result, cbp, in the whole range of 
fi. Of course, the agreement is not surprising for large 
/i, where replica symmetry holds. For smaller /i, the ob- 
served agreement indicates that replica symmetry break- 
ing corrections to the ground state energy are small (less 
than 1%). A similar conclusion was reached in Ref. l20l 
for the ± J distribution with zero mean. 

In particular, for jj. — Q we obtain Cqo — —1.38 ± 0.04 
(z = 4) and Coo = -1-72 ± 0.02 (z = 6), where the 
errors take into account the uncertainty on the correc- 



FIG. 2: Size dependence of the ground state energy, for z = 4 
and different values of fi. The lines represent the best fits 
with the form E/N = Cao + bN''^/^. The N = oo data (origin 
of the X— axis) are obtained in the BP approximation. 

tion exponent w, to be compared with our BP result 
EBP -1.351±0.002 (z = 4) and esp = -1.737±0.002 
(z = 6). It is also interesting to compare this with the 
ground state energy per spin found in twe^ and three 
dimensions'^^ (which have coordination number z = 4 
and z = 6, respectively) with Gaussian couplings and 
^ = 0, which is e^o = -1.31453(3) and e^o = -1.7003(1) 
respectively. 



B. Ground state magnetization 

In Figs. 01 and 0] the symbols show, for z = 4 and 
6 respectively, the average ground state magnetization 
m — [M]j, where M = J2i '^i ^^'^ [■■■].! denotes the 
sample average, as a function of fi for different system 
sizes A. The lines show the N = oo result in the BP 
approximation. For small /x, the magnetization vanishes 
as 1/Vn. For large /x, the finite- A^ data agree with the 
BP result within the error bars, with negligible finite-size 
corrections (again, we recall that the BP approximation 
is exact for sufficiently large /x, hence the agreement is 
expected). From the point at which the BP magnetiza- 
tion vanishes, we estimate the critical coupling strength 
Aic = 0.742±0.005 (z = 4) and fic = 0.546±0.005 (z = 6). 
Note that recursion relation Eq.lSJ admits two symmet- 
ric solutions for fi > /Xc- Hence, in the stochastic pro- 
cedure the magnetization will oscillate between positive 
and negative values, with an oscillation "time" (number 
of iterations) Mq that increases with the population size 
J\f and with /i. Therefore, to compute the magnetization 
correctly, we need A^o ^ To do this, we increased 
the size of the population progressively from A/" = 10'^ to 
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FIG. 3: Average ground state magnetization as a function 
of fi, for z = A. Symbols: branch-and-cut results (statisti- 
cal errors are smaller than the symbols). Line: BP results 
with a population size ranging from M = 10'^ (away from the 
transition) to M — 10^ (near the transition), and with and 
M — 10* iterations of the stochastic algorithm. 

N = 10^ as /i approached ^jLc- (Residual oscillations very 
close to /ic introduce a small systematic errop^Si, which is 
reflected in the errors for /ij. quoted above.) 

Another estimate of /ic can be obtained from the 
Binder cumulanljiS 

where [• • •] is now the "time" average. In the limit 
N oo, g(/i) = for /i < /ic and g(/i) = 1 for /i > /ic, 
hence (/(/i) can be used to locate /ic. As shown in Fig. [SJ 
the variation of the Binder cumulant with /i sharpens as 
M increases, an effect of the sign oscillations of the mag- 
netization, which become less important as A/" increases. 
From M = 10^ we estimate 

/if^ ^ 0.743 ± 0.005 [z^A) 
/if^ = 0.547 ± 0.005 (z = 6) 

which agrees with the above estimate from the average 
magnetization. We also verified that with these values of 
/ic, the magnetization obeys msp = a{ii — Hc)^ for /i ~ 
/ic, with the mean- field exponent (3 = 1/2 and a ~ 0.23. 

Klein et ali^ solved the BP recursion in the vicinity of 
/ic using the mean random field approximation (MRF). 
Their results /if = 0.775 (z = 4) and /if = 0.587 
{z = 6) (obtained after rescaling their value by an appro- 
priate normalization factor -y/i) are slightly larger than 
our result /i^^. 

In order to obtain an estimate of /ic from the finite- 
branch-and-cut data, we computed the Binder cumulant 



FIG. 4: Same as Fig. 01 but for 2 = 6. 
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FIG. 5: Binder cumulant from the stochastic solution of the 
BP ansatz, for three different sizes of the stochastic popula- 
tion TV. 

(/(/i. A''), defined as in Eg. 1)11(1 but with the time average 
replaced by the sample average. According to finite-size 
scaling, the curves for f/(/i, N) as a function of /i for var- 
ious A^ must cross at the critical point /i = /ic- In Fig.|^ 
we plot the Binder cumulant in the vicinity of the in- 
tersection point (note that the horizontal scale is much 
larger than that of Fig. [SJ, from which we obtain 

/ic = 0.77 ±0.02 (z = 4) 
/ic = 0.56 ±0.02 (z = 6). 

This agrees with within the error bars, suggesting 
that also for the magnetization replica symmetry break- 
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FIG. 6: Binder cumulant as a function of fi for various system 
sizes. Only the region around the phase transition is shown. 
The lines are only a guide to the eye. 
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FIG. 7: Scaling plot for the Binder cumulant. The symbols 
are the same as in the corresponding panels in Figure 15] Note 
the steeper shape of the scaling function for z — 6. 

ing corrections are small, causing a shift of fic of less than 
3 — 4%. Replica symmetry breaking corrections are ex- 
pected to increase with z. In the Sherrington-Kirkpatrick 
model (which is the z ~* oo limit of the present model), 
corrections shift from 1.25 to 1. Although our numeri- 
cal estimate of Hc is slightly larger tha.n instead, this 
could be a statistical fluctuation or a finite-size effect. 
The small size of replica-symmetry-breaking corrections 
to jjic suggests that the mixed ferromagnetic spin-glass 
phase is narrow for these values of z, as also recently 
indicated for the three-dimensional Ising spin glass^. 



FIG. 8: Scaling plot for the ground state magnetization. The 
symbols are the same as in Fig. |S| 

The Binder cumulant is expected to satisfy the follow- 
ing finite-size scaling relation^^ for /.t ~ /ici 

5(/i,iV)=5(^i/K^)(^_^j) (12) 

where d„ is the upper critical dimension, which for the 
Ising spin glass is c?„ = 6. As usual, by plotting g{^, N) 
against N^^^'^"'^\fj, — fic) with correct parameters fic and 
ly, the data points for different system sizes should col- 
lapse onto a single curve near (/i — fic) — 0. As shown 
in Fig. 13 using the estimates of obtained above and 
the mean- field exponent v = 1/2 we obtain a good data 
collapse, showing that finite size scaling is well satisfied 
in our range of sizes. 

In Fig. |S1 we also show scaling plots for the average 
magnetization m(/i, iV) = [i\/],7, whose scaling form is 

with the mean field exponent (3 = 1/2. The data show a 
good scahng collapse for ^ < fic- 

VI. TYPICAL RUNNING TIME OF OUR 
BRANCH AND CUT ALGORITHM 

In this section we study the running time of our pro- 
gram as a function of the mean coupling strength /i. In 
computer science the complexity of a problem is classi- 
fied in terms of the worst- case running time of its solution 
algorithmaSi^. Central notions here are the complexity 
classes P and NP. Informally, the class P consists of all 
decision problems (namely, problems whose solution can 
only be "yes" or "no" ) for which at least one algorithm is 
known that can generate an answer in polynomial time, 
even in the "worst case". The class NP consists of all 
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decision problems for which, if for a given instance the 
answer is "yes" , then there is a certificate from which 
the correctness of the answer can be verified in poly- 
nomial time. For example, the question "Given a spin 
glass instance, is there a spin configuration with energy 
less than or equal to E'o?" belongs to NP. If for a given 
instance the answer is positive, then there is a spin con- 
figuration with correct energy, and its correctness can be 
verified in polynomial time. Only the existence of such 
a certificate (spin configuration, in the above example) 
is required, not the ability to find it in polynomial time. 
The class NP contains P, but it might be larger (many 
believe it is larger, and answering the question whether 
P = NP is an important open problem). NP- complete 
problems are the "most difficult" in the class NP, in the 
sense that no polynomial algorithm is known for solving 
them, and if a polynomial algorithm could be found for 
one of them, this would imply that all of them arc poly- 
nomially solvablc^^. The classes P and NP are defined 
for decision problems, but similar ideas apply to combi- 
natorial optimization problems as well. Informally, an 
optimization problem is called NP-/iard, if it is at least 
as difficult as every NP-complete problem. In particular, 
an optimization problem is NP-hard if the associated de- 
cision problem is NP-complete. This is true for many 
optimization problems, e.g. the maximum cut problem 
or the travelling salesman problem. 

In practice, the running time can vary greatly from an 
instance of the problem to another, and the worst-case 
running time might very rarely occur. Recent work has 
therefore focused on the average running time with re- 
spect to random instances drawn from some probability 
distribution. Instead of the average one can also ana- 
lyze the median, or typical running time, which has the 
advantage of being less infiuenced by the occurrence of 
exponentially rare samples with huge running times. 

It should be noted that, unlike the worst-case complex- 
ity classification discussed above, which is an algorithm- 
independent feature of the problem itself, in general the 
typical running time can be different for different algo- 
rithms and implementations that solve the same problem. 

Returning to our problem, as mentioned in Section III, 
finding the ground state of the Bethe-lattice spin glass is 
an NP-hard problem. For all values of /i the possible 
realizations of the disorder are the same as for fi = 0. 
Hence, the algorithm has an exponential worst-case run- 
ning time even on instances deep in the ferromagnetic 
phase. However, for large /i highly frustrated realizations 
are very unlikely to appear, hence the typical running 
time will decrease as fj, increases. The question we ask 
here is whether, for large N, the running time undergoes 
a sharp transition as a function of /i and, if so, whether 
the transition coincides with the spin-glass/ferromagnet 
phase transition. 

One may use the CPU time as a measure of the running 
time. However, the CPU time is machine-dependent, 
hence it is not suitable when different computers are used. 
Furthermore, it is hard to separate out the infiuence of 




FIG. 9: Median running time (measured in number of solved 
linear problems) as a function of fi for different system sizes. 
Inset shows the same for 2: = 6 and A'' = 200, 140, 100, 50. 

size-dependent hardware effects on the CPU time (for ex- 
ample, small problems can be fully stored in the cache 
and therefore run faster). To avoid these problems, in the 
following we use the number of linear problems solved, 
nips, as a measure of the running timeii. 

In Fig. 1^1 we show the median running time so de- 
fined as a function of for ^ = 4, 6 and different sys- 
tem sizes. Clearly, ground states are calculated quickly 
in the ferromagnetic region, while in the spin-glass phase 
the running time increases dramatically (note the loga- 
rithmic scale on the vertical axis), and is approximately 
constant within the entire spin glass phase. The variation 
becomes more pronounced as N increases, suggesting a 
sharp discontinuity in the iV — s- 00 limit around fi w 0.8 
{z = 4) and n ~ 0.6 {z = 6), which is close to the spin- 
glass/ferromagnet transition point /ic determined in Sec- 
tion El 

As shown in Fig. 1101 deep in the ferromagnetic phase 
the data is consistent with a polynomial increase of the 
running times with N. For smaller values of /i, the curves 
are bending upwards, indicating that the running time 
increases faster than any polynomial. This is also the case 
for fi — 0.8 (z = 4) and for fi = 0.6 {z = 6, not shown). 
Hence from this data, it seems that the change in the 
typical running time occurs at a value of ^ larger than 
/ic, although it is difficult to locate a precise transition 
point. A mismatch between phase transition and change 
of the running time has been observed before, e.g. for a 
simple algorithm solving vertex cover™. 

We have fitted the data in Fig. ^] with a function 
of the form nips{N) ~ exp(6A^'^). For = 0, we find 
b = 0.026(9), c = 0.87(5) for z = 4, and 6 = 0.007(3) and 
c ~ 1.24(8) for z — 6, but the data exhibits in both cases 
a considerable scatter around the fitting region, prohibit- 
ing to conclude in a definite way that the typical running 
time is exponential. Nevertheless, the data strongly sug- 
gest so. 
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FIG. 10: Median running time as a function of N for different 
coupling strengths p in logarithmic scale. The straight lines 
represent power-laws c * A'^'' with — 0.699 (z = 4, /i = 1.2), 
C = 0.677 [z = 4,fi = 1.6) and C = 0.709 {z ^ 6,n = 1.6), 
respectively, showing that in the ferromagnetic phase the me- 
dian running time is polynomial. 

VII. CONCLUSIONS 

We have studied the ground state of a diluted mean- 
field Ising spin glass model with fixed connectivities 
z = 4, 6 and Gaussian distribution of the couplings, with 
mean fi and unit variance. We have apphed a branch- 
and-cut algorithm, a sophisticated technique originat- 
ing in combinatorial optimization which guarantees to 
find exact ground states. Our motivation was to study 
the spin-glass/ferromagnet transition and relate it to the 
change in the typical running time of our algorithm. 

From the study of the Binder cumulant, we have 
obtained values for the critical coupling strength, /ic. 
We have also solved the model in the Bethe-Peierls 
approximation, using an iterative stochastic procedure. 
In this approximation we obtain a critical coupling 
strength, fj,^^, which agrees with the branch- and-cut 
estimate within the error bars of the latter, indicating 
that replica symmetry breaking effects are quantitatively 



small. Finite-size scaling is well satisfied for systems of 
size larger than iV « 30. 

We have also analyzed the ground state energy, and 
shown that the branch-and-cut results, extrapolated 
to the thermodynamic limit, are in very good agree- 
ment with the Bethe-Peierls results, again indicating 
that replica symmetry breaking effects are quantitatively 
small. In the spin glass region, finite-size corrections are 
well described by a N~^^^ dependence. 

We have investigated the typical running time of our 
implementation of the branch-and-cut algorithm, which 
we defined as the median number of linear programs 
needed to find the ground state, with respect to a uni- 
form distribution over the space of instances. We have 
shown that finding ground states is "hard" in the spin- 
glass phase, and "easy" deep in the ferromagnetic region, 
with a sharp variation at a value of /i slightly larger than 
fic- The data indicate that while the worst-case running 
time is always exponential in the system size, the typical 
running-time is polynomial in the ferromagnetic phase 
and super-polynomial in the spin-glass phase. 

Our understanding of what makes a problem compu- 
tationally hard is still very weak. In this paper, we have 
shown that in a standard hard problem from physics, 
the Ising spin glass, a "physical" phase transition has a 
dramatic effect on the performance of a solution algo- 
rithm. Although in principle the "typical hardness" is 
algorithm-dependent, it is reasonable to expect that the 
phase transition will influence to some extent the running 
time of many other solution algorithms. Furthermore, we 
expect that similar phenomena occur in other well-known 
physical models. 
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